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Abstract 

In this paper, we study the exact solution of the restricted iso- 
chemical model of H2 molecule with fixed nuclei recently proposed 
by Santilli and Shillady in which the two electrons are assumed to 
be bonded/correlated into a quasiparticle called the isoelectronium. 
Under the conditions that: 1) the isoelectronium is stable; 2) the 
effective size of the isoelectronium is ignorable, in comparison to in- 
ternuclear distance; and 3) the two nuclei are at rest, the Santilli- 
Shillady model of the H2 molecule is reduced to a restricted three- 
body system essentially similar to a neutral version of the H2 ion, 
which, as such, admits exact solution. Our main result is that the 
restricted three-body Santilli-Shillady approach to H2 is capable to fit 
the experimental binding energy, at the isoelectronium mass param- 
eter M = 0.308381me, although under optimal internuclear distance 
about 19.6% bigger than the conventional experimental value, indi- 
cating an approximate character of the three-body model. 







1 Introduction 



In this paper, we study isochemical model of the H2 molecule recently in- 
troduced by R. M. Santilli and D. D. Shillady which is characterized 
by the conventional H2 model set up plus a short-range attractive Hulten 
potential interaction between the two electrons originating from the deep 
overlapping of their wave functions at mutual distances of the order of 1 fm; 
see also 0. If one assumes that this attractive potential is strong enough 
to overcome Coloumb repulsion between the two electrons, they can form 
electron-electron system called isoelectronium. The isoelectronium is char- 
acterized by "bare" mass M = 2me, as a sum of masses of two constituent 
electrons, charge — 2e, radius about 10~^^ cm, and null magnetic moment. 
The used Hulten potential contains two real parameters, one of which is the 
isoelectronium correlation length parameter rc, which can be treated as an 
effective radius of isoelectronium. 

The main structural difference between the Santilli-Shillady isochemical 
model and the conventional quantum chemical model of the H2 molecule, is 
that the former admits additional nonlinear, nonlocal, and nonpotential, thus 
nonunitary effects due to the deep overlapping of the wavepackets of valence 
electrons at short distances, which are responsible for the strong molecu- 
lar bond. In a first nonrelativistic approximation, Santilli and Shillady [|l| 
derived the following characteristics of the isoelectronium: total rest mass 
M = 2me, charge — 2e, magnetic moment zero, and radius 6.84323 x lO^^^cm. 
The value M = 2me of the rest mass was derived via the assumption of a 
contact, nonpotential interactions due to the mutual wave-overlapping suffi- 
ciently strong to overcome the repulsive Coulomb force. The nonpotential 
character of the bond was then responsible for the essential lack of binding 
energy in the isoelectronium, and the resulting value M = 2m f.. However, the 
authors stressed in that the isoelectronium is expected to have a non-null 
binding energy, and, therefore, a rest mass smaller than 2m f.. One argument 
presented in [|I| is that, when coupled in singlet at very short distances, the 
two electrons eventually experience very strong attractive forces of magnetic 
type, due to the two pairs of opposing magnetic polarities, resulting in a 
bond. The potential origin of the bond then implies the existence of a bind- 
ing energy, resulting in a rest mass of the isoelectronium smaller than 2me- 
Also, in the subsequent paper P, Santilli pointed out that the isoelectronium 
can at most admit a small instability. 
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As a result of a correlation/bonding between the two electrons, Santilli 
and Shillady were able to reach, for the first time, representations of the 
binding energy and other characteristics of H2 molecule which are accurate to 
the seventh digit, within the framework of numerical Hartree-Fock approach 
to H2 molecule viewed as a four-body system with fixed nuclei, and with the 
use of Gaussian screened Coloumb potential taken as an approximation to 
the Hulten potential 

On the other hand, the above mentioned strong short-range character of 
the electron-electron interaction suggests the use of approximation of sta- 
ble isoelectronium of ignorahly small size, in comparison to the internuclear 
distance [|l|. Indeed, under these two assumptions one can reduce the con- 
ventional four-body structure of the H2 molecule to a three-body system (the 
two electrons are viewed as a single point-like particle). Furthermore, in the 
Born-Oppenheimer approximation, i.e. at fixed nuclei, we have a restricted 
three-body system, the Shrodinger equation for which admits exact analytic 
solution. 

So, we have the original four-body Santilli-Shillady model of H2 molecule, 
and the three-body Santilli-Shillady model of H2, which is an approximation 
to it. The former is characterized by, in general, unstable isoelectronium and, 
thus, sensitivity to details of the electron-electron interaction, while the latter 
deals with a single point-like particle (stable isoelectronium of ignorable size) 
moving around two fixed nuclei. 

Clearly, the three-body Santilli-Shillady model of H2 molecule can be 
viewed as H2 ion like system. For the sake of brevity and to avoid confusion 
with the H2 ion itself, we denote H2 molecule, viewed as the restricted three- 
body system, as H2. Note that H2 is a neutral H2 ion like system. 

The quantum mechanical problem of the restricted H2 ion like systems, 
associated differential equation, and its exact analytic solution have been 



studied in the literature by various authors |q]-[p|. 

In this paper we present the exact analytic solution of the above indi- 
cated restricted three-body Santilli-Shillady isochemical model of the hydro- 
gen molecule, study its asymptotic behavior, and analyze the ground state 
energy, presenting numerical results in the form of tables and plots. Our 
analysis is based on the analytical results obtained for thoroughly studied 
H2 ion. 

In Sec. 1^, we review some features of the four-body Santilli-Shillady model 
of H2 necessary for our study, and introduce our separation of variables in 
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the Schrodinger equation under the assumption that the isoelectronium is a 
stable quasiparticle of ignorable size. 

In Sec. ^, we review the exact analytic solution of the H2 ion like systems 
(which includes the H2 system) , and study their asymptotic behavior at large 
and small distances between the two nuclei. 

In Sec. ^, we use the preceding solution to find the binding energy of 
H2 system. We then develop a scaling method and use Ritz's variational 
approach to check the results. Both the cases of the isoelectronium "bare" 
mass M = 2me and of variable mass parameter, M = rjirie have been studied. 
All the data and basic results of this Section have been collected in Table ^ 

In Sec. 1^, we introduce a prehminary study on the application of Ritz's 
variational approach to the general four-body Santilli-Shillady model of 
where the isoelectronium is an unstable composite particle, in which case 
the model re-acquires its four-body structure, yet preserves a strong bond- 
ing/correlation between the electrons. 

In the Appendix, we present the results of our numerical calculations of 
the ground state energy of H2 ion and of H2 system, for different values 
of the isoelectronium mass parameter M, based on their respective exact 
solutions, in the form of tables and plots. 

Our main result is that the restricted three-body Santilli-Shillady isochem- 
ical model of the hydrogen molecule does admit exact analytic solution capable 
of an essentially exact representation of the binding energy, although under 
internuclear distance about 19.6% bigger than the conventional experimental 
value. The mass parameter M of isoelectronium has been used here to fit the 
experimental value of the binding energy, with the result M = 0.308381me 
(i.e. about 7 times less than the "bare" mass M = 2me). In this paper, we 
assume that some defect of mass effect may have place leading to decrease 
of the "bare" mass M = 2me. 

We also note that the value M = 0.308381me implies a binding energy of 
about 1.7 MeV, which is admittedly rather large. Recent studies by Y. Rui 
[|] on the correct force law among spinning charges have indicated the exis- 
tence of a critical distance below which particles with the same charge attract 
each others. If confirmed, these studies imply that the repulsive Coulomb 
force itself between two electrons in singlet coupling can be attractive at a 
sufficiently small distance, thus eliminating the need to postulate an attrac- 
tive force sufficiently strong to overcome the repulsive Coulomb force. As a 
result, a binding energy in the isoelectronium structure of the order of 1.7 
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MeV cannot be excluded on grounds of our knowledge at this time. 

Clearly, however, that due to the current lack of dynamical description 
of the above mentioned defect of mass, and the obtained result that the 
predicted internuclear distance is about 19.6% bigger than the experimental 
value, our study is insufficient to conclude that the isoelectronium is perma- 
nently stable, and one needs for additional study on the four-body Santilli- 
Shillady isochemical model of H2, which is conducted in a subsequent paper 
by one of the authors 0. 



2 Santilli-Shillady model of H2 molecule 
2.1 General equation 

The Santilli-Shillady iso-Shrodinger's equation for H2 molecule with short- 
range attractive Hulten potential between the two electrons can be reduced 
to the following form pi: 
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where Vq and Tc are positive constants, and R is distance between nuclei a 
and h. By using vectors of center-of-mass system of electrons 1 and 2, fa and 
Tb, originated at nuclei a and 6, respectively, we have 
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(for electrons we have mi = m2 = me). The Lagrangian of the system can 
be written 



C 



^ + ^ - V{n2) - W{r,a, ri,, r2a, r^,, R), (2.4) 
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Here, V is the potential energy of interaction between the electrons 1 and 2, 

V(ri2) = Vo- ^, (2.5) 

and is the potential energy of interaction between electrons and nuclei, 
and between two nuclei, 

g2 g2 g2 g2 g2 

Vr(ria,rib,r2a,r26,i?) = + — . (2.6) 

ria rib R 

Notice that ri = ria = rif,, and r2 = ''^2a = ''^26, because ri^ = rif, + ^ and 
= "^26 + -R, where R is constant vector. Similarly, 

= + i?, = ■ , n = ■ . (2.7) 

mi + 772,2 mi + m2 



Then, Lagrangian (|2.4|) can be rewritten as £ = C{ra, r^, ri2), 

/: = ^ + ^ - V^(ri2) - l^(r„, n, ri2, i?). (2.8) 

Here, M = mi +m2 is the total mass of the electrons, and m = mim2/ (mi + 
m2) is the reduced mass. Corresponding generalized momenta take the form 

^ dC ^- ^ dC 

PM = -^ = MrA. P^m = -^ = mfu. (2.9) 
OVA ori2 

The system reveals axial symmetry, with the axis connecting two nuclei. 
Also, for identical nuclei we have reflection symmetry in respect to the plane 
perpendicular to the above axis and lying on equal distances from the two 
nuclei. 



2.2 Separation of variables 

Santilli and Shillady [|I|] then assume that, as a particular case under study 
in this paper (not to be confused with the general four-body case), the two 
valence electrons of the H2 molecule can form a stable quasi-particle of small 
size due to short-range attractive Hulten potential, such that 

ri2 < ri2 < n. (2.10) 
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Therefore, we can ignore in Eqs.( p.2| ) and ( p.3| ), 



ria ~ ~ Ta, Tib ^ ^ Tfe. (2.11) 

The Hamiltonian of the system then becomes 

H=^ + §1 + nri2) + W{ra, n, R), (2.12) 

where 

2p2 

W^(r,,n,i?) = -— - — + (2.13) 

ra rb R 

In this approximation, it is possible to separate the variables r^^b and ri2- 
Namely, inserting |0) = ip{ra,rb, R)x{i^i2) into the equation |jl| 

/ a' „2 ^'„2 e-^^'/'^ 2e2 26^ e^, ^, 

-^V^, - —V 2 - K)^ —7^ + +R <^ = ^ 

y 2M 2m 1 — e ^i^/rc y 

(2.14) 

we obtain 

- Itt^ - + Vin2) + Wir^, r„ i?) - E = 0. (2.15) 

By separating the variables, we have the following two equations: 

-^'^uX + Vin2)x = ex- (2.16) 

and ^ 

In this way, under approximation (|2.1CI|) , the original four-body problem is 
reduced to a three-body problem characterized by two differential equations: 

1) Equation (|2.16| ), which describes the electron-electron system forming 
the bound quasi-particle state called isoelectronium, with "bare" total mass 
M = 2me and charge — 2e. This equation will not be studied in this paper, 
since we assume that the isoelectronium is permanently stable. 

2) Equation ( p.l7| ), which is the structural equation of the restricted three- 
body Santilli-Shillady isochemical model H2, in which the stable isoelectro- 
nium with "bare" mass M = 2me, charge — 2e, null magnetic moment and 
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ignorable size orbits around the two nuclei, hereon assumed to have infinite 
mass (the Born-Oppenheimer approximation). 

This paper is devoted to the study of the exact analytic solution of the 
latter equation, and its capability to represent the experimental data on 
the binding energy, bond length, and other characteristics of the hydrogen 
molecule. 



3 Exact solution for H2 ion like system 

In this Section, we present analytical solution of the Schrodinger equation 
for H2 ion-like systems in Born-Oppenheimer approximation, we analyze 
the associated recurrence relations, and asymptotic behavior of the solutions 
at large and small distances between the two nuclei. As it was indicated 
[HI], this problem arises when Santilli-Shillady model of H2 is reduced to the 
restricted three-body problem characterized by Eq. ( p.l7| ), which possesses 
exact solution under appropriate separation of variables. 



3.1 Differential equations 

In Born-Oppenheimer approximation, i.e., at fixed nuclei, the equation for 
H2 ion-like system for a particle of mass M and charge q is 

V^iIj + 2M(E + — + -)iIj = 0. (3.1) 
ra n 

In spheroidal coordinates, 

X = — - — , 1 < X < 00, (3.2) 
R 

2/=^^^^^, -1<2/<1, (3.3) 
< < 27r, (3.4) 
where i? is a fixed separation distance between the nuclei a and 6, and 

V72 ^ ( ^ ( 2 \ ^ 'i\ ^ \ 

^ ^ R\x^-y^) \di^'' ~^^di^d^^^~y^d^) ^^-^^ 

1 

+ 



i?2(a;2 - 1)(1 -y2)5<^2- 
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We then have from Eq.(prT 



dx^"^ 'dx dy^ ^ ' dy A{x^ - - y^) dip^ ^ ' 



H {x^ - r) + 2MqRx 



ij = 0. 



Here, we have used 



1 1 4 X , 
- + -=o^ 2- 3.7 

Obviously, Equation ( p.6|) can be separated by the use of the representation 

^ = f{x)g{y)e^^^, (3.8) 

under which we have two second-order ordinary differential equations, 



d ( , . , d \ /, MER^ „ \ ^ , ^ 

- 1)—/ 1 - [a - 2MqRx ^x^ + j / = 0, (3.9) 



dx \ dx' 



d ( ,^ 0. d \ / MER^ 



m 



2 



where A is a separation constant (cf. [§]). So, the problem is to identify 
solutions for / and g. 

3.2 Recurrence relations 

By introducing the re-formulations 

f^{x'-ir"f, (3.11) 

g^{l-y^rl^g^ (3.12) 

to handle singularities at a; = ±1 and ?/ = ±1 in Eqs.(^.9D and ( |3.1CI| ), 
respectively, we reach the following final form of the equations to be solved: 

(x^ - l)f" + 2(m + l)f - (A + m(m + 1) - Sx - cV)/ = (3.13) 

and 

(1 - y'^)g" - 2(m + 1)^' + (A - m(m + 1) - cV)^ = 0, (3.14) 
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where we have denoted 

, MER^ 

= — - — , a = 2MqR. (3.15) 

We shall look for solutions in the form of power series. Substituting the 
power series 

/ = E/fc^'' (3-16) 
9 = T.9ky\ (3.17) 
into Eqs. ( p.l3|) and (|3.14| ), we obtain the recurrence relations, 



c'/n-2 + aU-1 -{X-{m + n){m + n + l))/„ (3.18) 
-(n+l)(n + 2)/„+2 = 

and 

c^9n~2 -{X-{m + n){m + n + 1))^„ -{n + l){n + 2)^„+2 = 0, (3.19) 

from which coefficients fk and Qk must be found. Here, /o and go are fixed by 
normalization of the general solution. Note that the recurrence relation ( p.l^ ) 
contains term 2MqRfn-i raised from the linear term 2MqRx in Eq.( p.9|) . 

In the next two Sections we consider some particular cases of interest 
prior to going into details of the general solution. These particular solutions 
are important for the study of the general case. 

3.3 The particular case R = 

In the particular case R = 0, the two nuclei are superimposed, so that the 
system is reduced to a helium-like system. 



From recurrence relations ( ^.18 ) and ( |3.19| ) we obtain the following particular 
recurrence sequences, 

(A - (m + n){m + n + l))/„ - {n + l){n + 2)U+2 = (3.22) 
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and 

(A - (m + n)(m + n + 1))^„ -{n + l){n + 2)c/„+2 = 0, (3.23) 

which are equivalent to each other, and can be stopped by putting the sepa- 
ration constant 

A = (m + n)(m + n + l) =/(/ + !), (3.24) 

with m = This gives us well known solution for g in terms of 
Legendre polynomials, 

9 = a-y'r^'^Piiy), (3.25) 



where m = \m\, and 

1 d 
¥11 Jy 

The solution is the well known spherical harmonic function 



Yim = NirnPr{y)e'"''', (3.27) 
with normalization constant 



N, 



Im 



(/-m)!(2/ + l) 
\\ (/ + m)!47r 



(3.28) 



This solution corresponds to the case of an ellipsoid degenerated into a sphere, 
and we can put y = cos 6 for identification with the angular spherical coordi- 
nates {6, if). Equation in x corresponds to the radial part of the well known 
solution expressed in terms of Laguerre polynomials. 

3.4 The particular case q = 

In the particular case of zero charge, g = 0, we have from Eqs.( |3.9| ) and 

(03) 



2 



lUl-y^)±,)^U-cY--^A,^^. (3.30) 
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One can see that these equations originate straightforwardly also from the 
standard wave equation V'^ip+k'^ip = 0, in the spheroidal coordinates (x, y, ip). 
Recurrence relations ( |3.18D and ( ^.191 ) then become 



c'/n^2 -{X-{m + n){m + n + !))/„ - (n + l)(n + 2)/„+2 = 



and 



c^gn~2 - {\ - {m + n){m + n + l))gn - (n + l)(n + 2)g^ 



n+2 



0, 



(3.31) 



(3.32) 



which are equivalent to each other. 

A general solution for / is given by linear combinations of radial spheroidal 
functions R^^{c,x) of first, p = 1, and second, p = 2, kind 0, 



g (2m + r)!^_ 



m/2 



T- = 0,1 



T : 



(3.33) 



X J2 ^ 

r=0,l 



r+m—n 



(2m + r) 



Tl 



m+r 



[cx), 



where, 



— yn+l/2{z), 



(3.34) 
(3.35) 



and Jn+i/2{z) and Vn+i/2(-2) are Bessel functions of first and second kind, 
respectively. The sum in ( p.33|) is made over either even or odd values of r 
depending on the parity of ri — m. Asymptotics of R^^l^{c,x) and R^\{c,x) 
are 

^ (3.36) 
X) - — y — sm cx — -(n+i)n . (3.37) 



(c 



1 








COS 


cx 


cx 






1 








sin 


cx 


cx 







2 
1 



Particularly, to have well defined limit at x = we should use only spheroidal 
function of first kind, i?^]j(c, x), because Bessel function of second kind, 
Yn{z), has logarithmic divergence at z = 0. 
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General solution for g is given by linear combination of angular spheroidal 
functions of first and second kind , 

oo' 

S^r^lic,y)= E (3.38) 

r=0,l 
oo' 

S^2ic,y)= E dric)QZ+riy), (3.39) 

r=— oo 

where P^iy) and Q^niv) the associated Legendre polynomials of first and 
second kind, respectively. 

Expressions for radial and angular spheroidal functions, and correspond- 
ing eigenvalues A, for particular values of m and n, are presented in Ref. 

i. 

Coefficients d'^"'{c) are calculated with the help of the following recurrence 
relation: 

akdk+2 + iPk - Amn)4 + lkdk-2 = 0, (3.40) 

where 

_ {2m + k + 2){2m + k + ly 

~ (2m + 2A; + 3)(2m + 2fc + 5)' ^' ' 

2(m + k)(m + k + 1) - 2m^ - 1 , 

3k = (m + k)(m + k + l)+ \ '\ ' — cl 3.42 

Hk \ ^ )\ ^ ^ )^ (2m + 2A;- l)(2m + 2A; + 3) ' ^ ^ 

k{k - ly 

(2m + 2A; - 3)(2m + 2A; - !)■ ^ ' ' 

The calculation is made by the following procedure. First, one calculates 
NT. 

NT^2=lT-^mn-^ (r>2), (3.44) 

r 

NT = lo - ^mn. NT = iT - A™„, (3.45) 

1 Am? — 1 

7-= (m + r)(m + r+l) + -c2 1 - , (^ > 0). 

2 [ (2m + 2r — l)(2m + 2r + 3) J 

(3.46) 

Second, one calculates the fractions do/d2r and di/d2p+i with the use of 

_c[o_ _ dQ^(h d2r-2 
^2r '^2 d^ C?2r 
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and 

j^rm^ {2m + r){2m + r - ly 

{2m + 2r -l){2m + 2r + l)dr-2 

The coefficients do, for even r, and di, for odd r, are determined via the 
normahzation of the solution. 



3.5 The general case 

In this Section, we consider the general solution of our basic equations (|3.9| ) 



and ( |3.1CI| ). To have more general set up, we consider the case of different 
charges of nuclei, Zi and Z2- This leads to appearance of additional linear in 
y term in Eg. ( ^-101 ), so that both the ordinary differential equations become 
of similar structure. Also, we restrict consideration by analyzing discrete 
spectrum, i.e. we assume that the energy E <Q. 
Let us denote 

p=:^yZ2E, a = R{Z2 + Z^), h = R{Z2-Z^). (3.50) 

Then, Eqs. (|3.9| ) and (|3.1CI|) , for the general case of different charges of nuclei, 
can be written as 

^((a:^-l)^/™.(p,a;x)) 
+ (-AS - v\x^ - 1) + ax - f^k{v, a; x) = 0, (3.51) 

+ (^X^y] - p'{l - y') + by- g^,{p, b; y) = 0, (3.52) 

where we assume that the solutions obey 

|/„fc(p,a; 1)1 < oo, lira fmk{p, a; x) = 0, \gmqip,b; ±1)\ < oo. (3.53) 
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The eigenvalues A in Eqs.( p.51"D and ( |3.52| ) should be equal to each other, 

A;:i(p,a) = AL^J(p,6). (3.54) 



The general solution ijj{x,y,(p) of Eq. (|3.6|) is represented in the following 
factorized form: 

'4'kgm{x, y, v?; R) = Nkgm{p, a, b)fmkip, CL] x)gmq{p, h] y) -= . (3.55) 



The normalization coefficients N^grnip, b) in Eq. (|3.55|) are represented with 
the help of derivatives of the eigenvalues, A^],(p, a) and X^l{p,b), namely. 



^kqm{Pi a,b) — -j^ 



d\%{p,b) dXtl{p,a 



mq \ 

1 



dp dp 



(3.56) 



For a given indices k, q, m, and fixed values of Zi, Z2, and R, the discrete 
energy spectrum E can be determined from Eq.( p.54D . This equation has 
unique solution, p = pkqm{a, b). Then, by solving the relation stemming from 
(» 

PkqURiZ2 + Zi), R{Z2 -Zi)) = (3.57) 

in respect to E, we can find the discrete spectrum of energy, 

^j{R) = EkqmiR, Zi, Z2). (3.58) 

Number of zeroes, k, q, and m, of the functions g{y), f{x), and exp ±imip are 
the angular, radial and azimuthal quantum numbers, respectively. However, 
instead of k, q, and m one can use their linear combinations, namely, = 
k + q + m + 1 is main quantum number and / = g + m is orbital quantum 
number. 

To construct the general solution u{z), which is called Coloumh spheroidal 
function ^ (csf), in terms of angular CSF giy) and radial CSF /(a;), let us, 
again, use the form which accounts for singularities at the points z = ±1 and 

z = 00, 

u{z) = (1 - z2)W2 exp[-p(l ± z)]viz). (3.59) 

Then, we represent v{z) as an expansion, 

00 

^(^) = ^s{p, b, \)ws{z), (3.60) 

s=0 
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in some set of basis functions Ws{z). 

Now, the complexity of the recurrence relations depends on the basis. In 
the preceding sections, where the particular cases, R = and q = 0, have 
been considered, we used a power series representation. One can try other 
forms of the representation as well. For a good choice of the basis functions 
Ws{z), we can obtain three-term recurrence relation of the form 



s-l 



0, 



(3.61) 



where a^, jSg, and 7^ are some polynomials in p, b, and A. Then, using the 
tridiagonal matrix A consisting of the coefficients a^, jSg, and 7^ entering 
Eq.( |3.61| ), we can write down the equation to find out eigenvalues A^].(p, a) 
and A^^(p, 6). Namely, 



det A = F{p, b, A) = 0. 



(3.62) 



The matrix A has a tridiagonal form. This leads directly to one-to-one cor- 
respondence between det A and the infinite chain fraction, 



F{p,b,X)=Po- 



ao7i 



Pi 



ai72 



P2 



Pn+i — ■ ■ 



Qo7i «i72 
Pi- P2- 



_ Qn 
Pn' 

(3.63) 

In numerical computations, this relation allows one to find out eigenvalues 
A in an easier way due to simpler algorithm provided by the chain fraction. 
Consequently, one can compute the energy and coefficients of the expan- 
sion of eigenfunctions g{y) and f{x) by using the chain fraction. 

The result of this approach in constructing of the solutions depends on the 
convergence of the chain fraction. Analysis of the convergence can be made 
from a general point of view. Sufficient conditions of the convergence of the 
chain ( p.63| ), and of the expansion ( |3.60| ), are the following two relations: 



Ols-lls 



Ps-lPs 



< 



O-s+l 



Ps_ 

2a, 



PI 



1/2 



(3.64) 



Further analysis of the convergence depends on specific choice of the basis 
functions Us{z). 
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(i) Series expansion, Vs{z) = z" . In this case, the radius, Z^,, of conver- 
gence is 

as 



hm 



(3.65) 



Particularly, when a^+i/a^ ^ at s ^ oo the series ( p.60| ) converges at any 



(ii) For the choice of basis function Vs{z) in the form of orthogonal poly- 
nomials, the sufficient condition for convergence of Fourier series ( ^.6U| ) is 



< 1 - -. 

s 



(3.66) 



Below, we consider separately angular and radial CSF entering the general 
solution. 



3.5.1 The angular Coloumb spheroidal function 

For the angular Coloumb spheroidal function (acsf), it is natural to choose 
the basis functions Vs{y) in the form of associated Legendre polynomials, 
P^miy)- Indeed, they form complete system in the region y G [—1,1], and 
reproduce ACSF at p = 6 = (see Sec. |3]^). Inserting of the expansion 

9m,iP,b;y) = J2^sPr+miy) (3.67) 

s=0 



into Eg. ( p.52[ ) entails five-term recurrence relation. However, this relation, 



which is sometimes used, is not so suitable as the three-term relation. This 
is because the determinant of the corresponding pentadiagonal matrix can 
not be represented as a chain fraction. Nevertheless, in the case 6 = 0, i.e. 
for Zi = Z2, this five-terms recurrence relation is reduced to two three-terms 
recurrence relations, separately for even (c_2 = 0, Cq = 1) and odd (c_i = 0, 
Ci = 1) solutions of Eq. (|3.52|) presented in previous Section. 

For the general case 6 7^ 0, the expansions of g{p, b; y), handling singular- 
ities at the points y = ±1 and y = 00, respectively, as considered by Baber 
and Hasse 0, are 

00 

9m,{p, b; y) = exp[-p(l + y)] ^ CsPT+miv). (3-68) 
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9n.,{p. h- y) = exp[-p(l - y)] c'^Pr+miv), (3-69) 

s=0 

These expansions yield three-terms recurrence relation, 

psCs+l - KsCs + 6sCs-l = 0, c_i = 0, (3.70) 
where the coefficients for the case of expansion (|3.68|) have the following form: 

{s + 2m + 1) [b - 2p{s + m + 1)] 



Ps 



2(s + m) + 3 

Ks = {s + m){s + m + 1) — X, 

^ s[b + 2p{s + m)] 
' ~ 2(s + m)-l 



(3.71) 



To estimate convergence of these expansions, one can use the above made 
estimation of the convergence, with the following replacements: — > ps, 
Ps ^ i^s, Is ^ Ss, and as ^ Cs- For the expansion (p.79|) we have 



Ps~l5s 



l^s—ll^s 



(3.72) 



i.e., at p > 1, convergence takes place only at s > 2p. We should take into 
account this condition when choosing minimal number of terms in the chain 
fraction ( |3.73| ) which is sufficient to calculate A, to a required accuracy. 

The recurrence relation for the coefficients c'^ of the expansion ( p.69| ) 
differs from that of Eq.( |3.70| ) by the replacement p —>■ —p in formulas ( |3.71| ). 
Clearly, this replacement does not change the form of the chain fraction. 



F^y\p,b,x) = Ko- 



PqSi pi52 

Kl- K2- 



(3.73) 



So, in both the cases, (|3.68|) and (|3.69|) , the eigenvalues A can be found from 
one and the same equation. 



F^y\p,b,X) = 0. 



(3.74) 



In practical calculations with the help of this algorithm^ the infinite chain 
fraction ( p.63| ) is, of course, replaced by the finite one, F^l-^{p, b, A), in which 
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one retains a sufficiently big number of terms. Typically, iV > 10 provides 
very good accuracy. So, the eigenvalues are computed as the roots of the 
polynomial Qn+i{p, b, A) of degree N + 1, namely. 

Such a representation allows one to exclude singularities, associated to zeroes 
of the polynomial P/v+i(p, b, A), from Eq. (|3.74|) . Further, from the definitions 
( |3.63|) and ( |3.75| ) we obtain the following recurrence relation for the polyno- 
mial Qkip, b, A): 

Qk+l = Qkl^N-k — Qk-lPN-k^N-k+l, Q-1 = 0, Qo = 1, (3.76) 

with the use of which one can find Qn+i- Here, the coefficients Rg, Ps, and Sg 
differ from that of Eq.( |3.71| ) by the factor (1 + k^)"^/^. This factor does not 
change the recurrence relation ( |3.7CI|) . However, it makes possible to avoid 
accumulating of big numbers at intermediate computations. Indeed, from 
Eq.( 3.71|) for Kg it follows that the leading coefficients of the polynomials 



Qk would behave as k'^'', for example, for A; = 4 we would have 4^^, if we 
would not made the above mentioned renormalization of the coefficients ps, 
and 5s- The eigenvalue is found as an appropriate root of the polynomial 
Qn+i{p, b, X^y^). Clearly, for big A^, there is no way to represent in general the 
roots of Qn+1 analytically so one is forced to use numerical computations. 

In the numerical computations, to pick up the appropriate eigenvalue 
A^^(p, b) among + 1 roots of the polynomial Qn+i it is necessary to choose 
some starting value of A. For example, one can put the starting value at the 
point p = b = 0, where A^^(0, 0) = {q + m){q + m + 1). The ffist step is to 
increase discretely p ^ p + Ap and 6 — >■ 6 + Ab beginning from the starting 
point p = b = 0, at fixed values of m and q, and the second step is to find 
^mqiP + b + Ab) with the help of Eq. (|3.74|) . Repeating these steps one 



mq \ 

can find A^^ numerically as a function of p and b in some interval of interest. 

Also, asymptotics of A which will be studied in Sec. |3.6| are of much help 
here to choose the appropriate root. For example, for 6 = and A^ = 5 we 
obtain numerically from the determinant of the tridiagonal matrix consist- 
ing of the coefficients defined by Eq. ( |3.7lD , with — on the main digonal, 
and Ps and 6s on the upper and lower adjacent diagonals respectively, the 
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polynomial, 



deti = 0.003A^-0.2A^+(0.2p2+5.5)A^-(6.3p^+56)A^ + (1.5/+66p2+231)A^ 

(3.77) 

-(19p^ + 226p2 + 277)A + / + 4V + 186/. 
Only one of its six roots has asymptotics, 



A|p^o = 0.667p2 - 0.0148p^ + O(p^), 



(3.78) 



which reproduces, to a good accuracy, the asymptotics (|3.124|) . So, this is 
the desired root to be used in subsequent calculations. Also, observe the 
decrease of the numerical coefficients at higher degrees of A which control 
the convergence. 

Note that, at p ^ 1, the ACSF is concentrated around the points y = ±1 
so that expansion ( 3.68|) converges slowly. In this case one uses another, more 
appropriate, expansions. 



g^,{p, h- y) = {l- vT'^ exp[-p(l + y)] ^ c.(l + yf 



s=0 



9^qiP,b;y) 



y 



2\m/2 



exp[-p{l-y)]J2c'sil-yy- 



(3.79) 



(3.80) 



s=0 



Evidently, expansion ( |3.79| ) converges faster in the region [—1,0] while the 
expansion ( p.80|) converges faster in the region [0, 1]. Here, the coefficients Cs 
of the expansion ( ^.79| ) obey the three-term recurrence sequence ( ^.70| ), with 

Ps = 2(s + l)(s + m + l), 

Ks = s{s + 1) + (2s + m + l)(2p + m) + 6 - A, (3.81) 
6s = b + 2p{s + m). 
It is remarkable to note that expansions ( p. 791 ) and ( |3.80| ) converge at any y, 



and the corresponding chain fractions ( |3.63|) converge at any p since 



2p 



Ps-l6s 



Ap 



(3.82) 



Similarly, the coefficients c'^ obey the same relation, with the replacement 
6 ^ -6 in Eq.(p:8T]). 
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In practical calculations, one can use a combination of expansions 
and ( 3.79|) . Namely, the procedure is: from expansion ( 3.68| ) one finds eigen- 



values while the eigenfunctions are calculated from to Eg. ( ^.791 ). Of course, 
both solutions (|3.79|) and (|3.8CI|) should be sewed, for example, at the point 



y = 0, because the recurrence relations do not determine, in this case, a 
general normalization of the coefficients Cg and c'^. Particularly, the sewing 
condition, which defines the normalization of Cg and c^, has the form 

Ec. = E4 (3.83) 

s=0 s=0 

To derive the asymptotics of ACSF and its eigenvalues we can use an 
expansion in Laguerre polynomials, 

oo 

g^,{p,b;y) = (1 - y^r/^ exp[-p(l ± y)] ^ c,Lr+^(2p(l ± y)), (3.84) 

s=0 



L^(z) = ^^^(e-^^"+-). (3.85) 



The insertion of this expansion into Eq. p.52| ) and the use of the differential 
equation for Laguerre polynomials, 

^^L-(^) + (1 - ^ + m)-^^K{z) + nL^{z) = (3.86) 

yield recurrence relation ( |3.70| ). For the case of positive sign in Eq. ( |3.84| ), we 
should put 

Ps = -{s + m + 1) { s + 1 + ^ 



5 



Ks = -(2s + m + l) ls + m + l + — -2p]+(s + m)(m + l) + b-X, (3.87) 
V J 

Sg = -s (s + m + ^y 
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3.5.2 The radial Coloumb spheroidal function 

The radial Coloumb spheroidal function (rcsf) obviously should be written 
in a form suitable to handle singularities at the points x = 1 and x = oo, 
namely, 

Uk{p, a; x) = {l- x^r^^ exp[-p{x - (3.88) 
So, the equation for /(x) takes the form 

{x'^-l)f"{x) + [-2p{x'^-l)+2{m+l)x]f'{x) + [-\+m{m+l)+2pax]f{x) = 0, 

(3.89) 

where we have denoted a = ^ — (m + 1). In the case when the expansion 

fix) = Y.asUs{x) (3.90) 

implies a three-terms recurrence relation, the eigenvalues A^],(p, a) can be 
found from the chain fraction equation, 

F(^)(j9,a;A) = 0. (3.91) 

Also, the expansion which is of practical use has been considered by Jaffe 
[p^ . In this case, the expansion series (|3.90|) becomes 

fix) = ix + irY.^sX', (3.92) 

where x = (x — l)/(x + 1) is Jaffe's variable. By inserting ( |3.92|) into the 
equation for the function f{x), we get recurrence relation ( p. 611 ), where the 
coefficients are 

as = {s + l){s + m + 1), 
f3s = 2s^ + {2s + m + l){2p - a) - a - m{m + 1) + A = (3.93) 
= 2s{s + 2p-a) - {m + a){m + 1) - 2pa + A, 
7s = (s — 1 — cr)(s — m — 1 — (t). 
Also, for Jaffe series expansion, we have 



Ps-lf3s 
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i.e., the chain fraction converges at p > 0. One can see also that the Jaffe 
expansion converges at any x. 

In addition, the function f{x) can be expanded in associated Laguerre 
polynomials. 



/(a;) = (x + ir^a,L; 

s=0 



+m 



X) 



X = 2p{x — 1). 



(3.95) 



In this case, the recurrence relation is of three-terms form, and the coefficients 
are 



a. 



— (s + m + 1) 



2p 



= -(2s + m+l) 



a 

2p 



[s + m + l){s — m — a) 



(s + m + 1) 



(3.96) 



+2p(2s + m + 1) - (s + m)(m + 1) - a + A, 



7. 



— s 



2p 



— [s + m) 



sis 



a). 



As to numerical computation of the eigenvalues X^^\ expansions (|3.92| ) and 
( p.95| ) are equivalent because the chain fraction depends, in fact, only on /^^ 
and as'js+i- Indeed, by comparing Eq. (|3.93|) and Eq. (|3.96| ), one can easily 
see that in both cases Ps and ag^js+i are the same. Evidently, it then follows 
that the associated chain fractions are equivalent to each other. 

However, we should note that Jaffe's recurrence sequence, in general, is 
more stable, while Laguerre expansion (|3.95|) is more suitable to find out the 
asymptotics of fmk{p, a; x). 

Also, we note that the associated "radial" polynomials Qn+i, the root 
A''^-' (p, a) of which should be found, contain a much bigger number of terms, in 
comparison to the "angular" case. So, practically finding of radial eigenvalues 
is much harder than that of angular eigenvalues. 

At equal charges of nuclei, Zi = Z2, the equation for g, and the recurrence 
relation for gk, are the same as they are in the particular case q = considered 
in Sec. WA - A general solution for g is then given by ACSF (|]38D and (|]3|), 
with coefficients c?™" given by recurrence relation ( |3.4CI| ). 

In the reminder of this Section we would like to note that, in general, 
solving the recurrence relations can be made equivalent to solving associated 
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ordinary differential equations by making the z transform. In many cases 
the z transform helps to solve recurrence relations. Namely, one defines the 
function 



Ziz) 



EQjg 

s=0 ^ 



(3.97) 



associated to the coefficients entering Eq. (|3.61|) viewed as a function of 
discrete variable s. For as, (3.s, and 7^ given by Eq.( p.93D we obtain from 
Eq- CT) 

;i - m)z'^ + 2{2p - cT - 1)2 + 2(T + m + ll Z'+ (3.98) 



z(z - i fZ" + 



+ (a + m)(m + 1) + 2pa - A 
For the coefficients Cs, we define 

Y(z] 



[m + a)a 



°° c 

s=0 ^ 



and for ps, Ks, and Sg given by Eq.( p.81| ) we obtain from Eq. ( |3.7CI| ) 
z{z - ifY" + \{t - 2)z^ + 2{2p -t + 1)z + t\ Y'+ 



+ 



[m(m + 1) — mT)z — 2p{m + 1) + (r + m)(m + 1) + 6 — A 



where we have denoted 



'm{rn — r + 1)cqZ = 0, 



T = \- m + I. 

2p 



(3.99) 

(3.100) 

Z- 

(3.101) 



If one has solved these differential equations for Z{z) and Y{z), then, by 
making the inverse z transform, one can find the expansion coefficients 
and Cg (and thus the general solution of the problem). 



3.6 Asymptotics of CSF and their eigenvalues 

To analyze the exact solution, which is of rather complicated nonclosed form 
(infinite series) given in the previous Sections, it is much instructive to derive 
its asymptotics, which can be represented in a closed form. In this Section, 
we present the asymptotics at large {R — >■ 00) and small {R — > 0) distances 
between the nuclei, with a particular attention paid to the ground state. 
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3.6.1 Asymptotics at /? ^ cxd 

For increasing distances R between the nuclei, at fixed quantum numbers k, 
q, and m, we have increasing values of the parameters p, a, and b, 

p^{-2Ef/'^R/2^oo, {Z2 + Zi)R^ oo, b ^ {Z2 - Zi)R ±00. 

(3.102) 

Let us introduce the notation 



Z2 + Z1 



2p v^^ZB ' 
and assume that o; ~ 1 and /? ~ 1. 



Z2 — Z\ 



2p ^/^2E 



(3.103) 



ACSF at i? ^ 00. 

Let us consider the asymptotic expansion of ACSF. In this case, the equation 
for the Whittaker function, M^,f^{y), builds ansatz around the poles y = ±1. 
Here, the solution is constructed in two overlapping intervals, P_ = [— l,|/i] 
and P+ = [1/2, 1], with y2 < yi- Then, the asymptotics of ACSF gmq{p, 2p/3; y) 
in the interval have the form 



9mq{p,'2pp;y) 



X 



r(m + i) 

M,,m/2 (2p(l + y) + 2{k + 13) In i^) 



1/2 



while in the interval it is 
9mq(p,'^pP]y) = 



2r («+^) 

[i + 0{p-% yev^ 



(3.104) 



2V[k' + 



l+m 



X 



r(m + i) 

M,,,m/i (2p(l - y) + 2{k' + /3) In i±a 



r [k' + if^) 



1/2 



X 



J [l + 0(p-^)], yeP+. (3.105) 

Here, the coefficients d- and d+ (d?. + = 1) are defined by the relations 

1/2 

(3.106) 



sin7r(2K' — m — 1^ 



sin7r(2fi; — m — 1) + sin7r(2fi;' — m — 1) 



X 
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xsgn 



— cos7r(K — (m + l)/2) 
sin7r(/t' — (m + l)/2) 

sin 7r(2K — m — 1) 



11 



1/2 



(3.107) 



sin7r(2K — m — 1) + sin7r(2K' — m 
RCSF at i? ^ oo. 

Now, let us consider the asymptotic expansion of RCSF. The replacements 
X —y, p — > —p, a —(3 convert the radial equation around the point x = 
1 to the angular equation around the point y = —1. Thus, the corresponding 
asymptotics of RCSF are directly related to the above found asymptotics of 
ACSF. 

The RCSF, normalized to the first order in p, has the form 



fmkip,'2pa;x) = — 
ml 



X 



2{k + +my. 

kUx^ - 



nl/2 



X 



i.,m/2 (2p(x - 1)) [1 + Oip-')]. (3.108) 

Since the first index of Whittaker function in Eg. ( p. 10^ ) is k = A; + (m+ 1)/2, 
the function can be expressed in terms of Laguerre polynomials. 



Energy at R ^ oo. 
In the limit R oo, the Coloumb two-center problem is evidently reduced to 
two separate problems of Coloumb centers, with the charges Zi and Z2. Each 
of the atoms, eZi and eZ2, is characterized by a set of parabolic quantum 
numbers, [n, rii, n2, m] and [n', n[, n'g, m], which are related to each other by 
the relations 

n = rii + n2 + m + 1, n' = n\ + n'2 + m + 1. (3.109) 

The number k of zeroes of RCSF coincides with the number ni, for the angular 
functions of the left center, eZi, and with the number n'^, for the angular 
functions of the right center, eZ2- 

A series expansion in inverse power of R can be obtained in the form 

E,.„.UZM ^ -|^-| + |||-^(6A=-„=«)+ (3.110) 
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Z rc' 

:[ZiA(109A2 - 39^2 - Qm^ + 59) - Z2n{XJn^ - SA^ - 9m'^ + 19)] + 



IGR^Zf ' 



where A = rii — n2, and £5 e are defined via the expressions, 



n 



3 



£5 = ^[n^(1065A^ - 594n2A2 + 1230A2 - 234m2A2 + 9m^+ (3.111) 

+33n^ - 18nW - ISm^ + 105 - 138^^) + 4n|A(21 A^ - nin^ + QSm'^ - 189)]. 



4 



71/ 

Eg = -^^[nzA(-2727A^ + 207671"^ A"^ - 5544A2 + 1056m2A2 - 93m^- 

^ (3.112) 
-273n^ + 78n^m'^ + A50m^ - 1533 + 147071^) + 2n|(-207A^ + 1044n2A2+ 

+2436A2 - 576A2m^ - 4272^ + 371 - 162m^ + 42mV - 89n^ + 15m^) + 

+2r2|A(3A2 - 69n^ - 117 - 33771^)], 

where 77^ = nZijZi. 

Eq. ( 13.1101 ) gives the multipole expansion in the electrostatic energy of the 
interaction between the atom eZi and the far-distant charge Z2 (so called 
eZi-terms). 

Note that expansion ( p. HQ ) can be obtained by ordinary perturbation 
techniques as well. Indeed, the degrees of Z\ display the orders of the mul- 
tipole moment of the atom eZ\. 

The series of terms corresponding to the other atom, 6^2, is obtained 
from Eq.( p.llOl ) with the use of self-evident replacements, Z\ <r^ Z2, n — > n', 
A A', and n2 — ^ 773. 

Finally, the energy of the ground state Isag of the molecular ion, for which 
Zi = Z2 = I (equal charges of nuclei), can be written, to a high accuracy, as 

1 9 15 213 7755 1733 86049 / 1 
1000 ( , - "2"4^~2^"4^"6lR8~ 2^~TaRio ~ V^ii, 

(3.113) 
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3.6.2 Asymptotics at i? — > 
Energy at i? — > 0. 

In the case of positive total charge, Z = Zi + Z2 > 0, and at i? ^ 0, we can 
use perturbative approach to ZieZ2 problem, without using a separation of 
variables. Namely, the Hamiltonian of the system Z\eZ2 is represented as 
the sum 

tl = Il^^^W = — -^~^. (3.114) 

2m ri r2 

The operator is usually chosen as the Hamilton operator of the so called 

united atom, 

which is placed on the z-axis at the point z = zo, 

The point (0, 0, Zq) is called center of charges due to the fact that it lies at 
the distances 

Ri = ^R and R2 = ^R, (3.117) 

from the left and right atoms, respectively. 

We choose a spherical coordinate system, (r^, -(^c, v^), with the origin at 
point (0, 0, 2:0), and the angle measured from z-axis. Then, the eigenstates 
iHnm of the operator H^^ are 

ili,t{rc) = RNi{rc)Yr{A,^), (3.118) 

while the eigenvalues are given by 

ENt = -^- (3.119) 

The matrix W^[^' of the perturbation operator W is diagonal on the func- 
tions ip^imi^c) of the atom if zq is defined by Eq.( |3.116| ). Below, the first two 
terms of the expansion of energy in powers of R are given, 

ENlm{Zi, Z2, R) = 
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2ZiZ2\l{l 



3m^ 



For the ground state of the ZieZi system, with equal values of the charges, 
one can find the following expression for the energy, up to the second order 
of perturbation: 



1 + -iZRf - -iZRf + -^(ZRf - —iZRflnZR + 

2 6^^ 6^ ^ 2160^ ' 36^ ' 



CSF at i? ^ 0. 
For fixed quantum numbers we have, at i? — > 0, 



(3.121) 



p = (-2^)1/2^/2 ^ 0, a = (Z2 + ^ 0, 

Let us denote 

2 Z2 + Zi 



b={Z2- Zi)R 0. 

(3.122) 



a 



2p y/^2E 



(T + m + 1, P 



2p V^2E 



(3.123) 



In this notation, the energy is ii^ = — Z^/(2q;)^. Let us consider asymptotics 
of CSF of the ground state of the molecular ion. 

ACSF at i? 0. 

The power series expansion of ACSF goo{p, 2pj3] y) in small parameter p can 
be obtained by expanding it in the Legendre polynomials. For the eigenvalue 
Xooip, 2p(3), we then get 



xl'J{p,2p(3) = il-0' 



lp'-,^p\i + np')+oip') 

3 135 



(3.124) 



RCSF at i? ^ 0. 

To expand RCSF /oo(p, 2p(l + <7);x), at p — > 0, a = O(p^), we use Jaffe's 
expansion. 



/oo(p,2p(l + a);x) = exp (-px)(l + x)'" ^ a,x', X 



X — 1 



(3.125) 



.=0 2; + 1 ■ 

where a^'s obey three-term recurrence relation with the coefficients ( 3.93| ). 
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For the eigenvalue X^oq{p, 2p{l + a)), we get 

aS^(p, 2p{l + a)) = a{l + 2p) + a\l + 4p\nAp-f) + o{p^). (3.126) 

RCSF of the ground state of the molecular ion, Zi = Z2 = I, can be presented 
as 

/oo(p, 2p{l + a); x) = exp {-px){l + x)'^[l + a'U,{x) + o{p% (3.127) 
where x = {^~^) / (x+l) is Jaffe's variable and Li2(x) is dilogarithm function, 



".(x) = ES = -/i"(i-«)f. 

n=l n ^ 



(3.128) 



The ground state energy is defined as a function of three parameters, Zi, Z2, 
and R, 

xlSip, 2pf3) = aS(p, 2p(l + a)) (3.129) 

Combining Eqs.( p.l2^ ), ( |3.126| ) and ( p.l29| ), we get series expansion for the 
ground state energy of the Z^eZ^ system in the form 



-E'ooo(^i5 ^2, -R) — 
-Z^ + '^-Z^Z^^ZRf - '^ZiZ^iZRf + '^Z,Z2 ( 1 



64^1^2 



{zrY- 



45 



27^2 

(ZR)^ + o{{ZR)^). (3.130) 



Comparing Eqs.( |3.121| ) and ( |3.130| ) we see that the terms proportional to 
{ZRY and (ZR)^ coincide. The next order corrections makes a difference; 
Eq. ( |3.121| ) obtained by the second-order perturbation is of less accuracy. A 
practically achieved accuracy of the first-order perturbation (|3.120|) and of 
Eq.( |3.130D is the same; at ZR < 0.1, the discrepancy is not bigger than 1%, 
and becomes sharply smaller with the increase of the parameter ZR. 



3.6.3 Quasiclassical asymptotics 

At -R ^ 00, for eZi solutions we have 

A^|,(p, 2pa) = —2p{2K — a) — k{2k — a — m) + 



(3.131) 
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3Ka + — m^] 



-2\ 



X'^alip, 2pf3) = 2p{2x + /3) - x(2x + 2/? - m)- 
2p 

Ax) 



(3.132) 



From the equahty A^^ = A^^, we get the expansion for Ej{R) which coincides 
with the asymptotics ( p.llO| ), up to terms of the order of R~^. 
In the hmit R ^ 0, the following expansions are justified, 



ink 



2p 



- (A; + 1/2) 



AS = (/ + l/2f + ^ 



+ 0{p' 



mq 



and we get, by using the equation A*^"^^ 



^(/ + 1/2)2 
A(^), 



0{p' 



ENlm{R) 



Zi 



N 



ZlZ2{Zi 



z,? 



AN^{1 + 1/2)5 



(3.133) 
(3.134) 



[(/ + 1/2)2 -3m2]. (3135) 



Expression ( p.l35| ) for the energy coincides with the asymptotics (|3.120| ), 
up to 0(/~2). Clearly, an accuracy of the quasiclassical Eqs. (|3. 131|) -( |3. 13^ ) 
becomes higher for a greater number of zeroes of the solutions, k and q. 
However, even for the ground state, k = q = m = 0, these equations give 
a good approximation for the energy in both limiting cases, R —>■ and 
R —y oo. Also, we note that corresponding numerical calculations showed 
that for the intermediate values of R the terms Ej[R) can be determined 
within the quasiclassical approach with accuracy of about 5%, or more |TT 



4 Scaling method and binding energy of three- 
body Santilli-Shillady isochemical model H2 

To find the ground state energies of H2 and H2, we use computations of 
the Isa terms of H2 ion and of H2 based on the above presented exact CSF 
solution by solving the corresponding equations A'-^-' = A*^^-*. 
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The angular and radial eigenvalues A are found as solutions of the equa- 
tions containing infinite chain fractions presented in previous Sections which 
should be interrupted and then solved numerically, to a required accuracy. 

Our primary interest is the study of Santilli-Shillady model H2 system. 
However, we present the results for H2 ion as well to check our calculations 
and to use them in the scaling method described below. 

Also, the reader should keep in mind that we are primarily interested in 
ascertaining whether there exist a non-zero value of R for which a fully stable 
and point-like isoelectronium permits an exact representation of the binding 
energy of H2 molecule. 

It should be noted that the ground state electronic energy is obtained 
as a function of the parameter R due to Eq.( |3.58D . By adding to it the 
internuclear potential energy 1/-R, we obtain the tota/ ground state energy of 
the system, so that at some value R = Rmin, the total energy E necessarily 
has a minimum, if the system is stable. This is the way to determine uniquely 
the internuclear distance under an exact representation of the total energy 
E. 

To have an independent check of the result for the total ground state en- 
ergy of H2 (with the stable and point-like isoelectronium) obtaining from the 
exact solution, we develop a scaling method based on the original Schrodinger 
equation for H2 ion like system. Namely, it appears that one is able to cal- 
culate the ground state energy as a function of R for any H2 ion like system 
with equal charges of nuclei, Zi = Z2, provided that one knows the ground 
state energy as a function of R for the H2 ion itself. It should be pointed out 
that the scaling method does not depend on the obtained solution because 
it reflects, in fact, the scaling properties of the Schrodinger equation itself. 

In addition, we use below Ritz's variational approach to H2 like systems 
to find out the approximate value of the ground state energy of H2, as well 
as to check the result provided by the exact solution, and to demonstrate the 
accuracy of the variational approximation. 

Our general remark is that in both approaches, the exact solution and 
Ritz's variational solution, we use Born-Oppenheimer approximation (fixed 
nuclei). Clearly, taking into account the first order correction, i.e., zero 
harmonic oscillations of the nuclei in H2 around their equilibrium positions, 
we achieve greater accuracy. 

But we still have a significant inaccuracy in the value of dissociation 
energy due to the fact that H2 system has the lightest possible nuclei (two 
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single protons). 

To estimate this inaccuracy, one can invoke Morse's potential customarily 
used for diatomic molecules. In particular, the analysis for H2 molecule shows 
that the ground state energy of harmonic oscillations of the nuclei receives 
1.4% correction due to the first anharmonic term. 



4.1 Exact representation of binding energies of H2 ion 
and three- body H2 system 

The CSF based computations for H2 ion were presented, for example, by 
Teller |]T^, Bates et al. |T^, and Wind |jT^, and we do not repeat this study 
here for brevity, while we shall just describe the procedure and present our 
final numerical results in Appendix, Table ^. 

In particular. Teller presented a plot of the resulting function Eisa{R) 
which provides a good accuracy, and Wind used the exact solution to present 
a table of energy values in seven decimal places for distance values of R up to 
20 a.u. in steps of 0.05 a.u. However, these results cannot be used directly in 
our case since the repulsive potential between the nuclei, 1/R, has not been 
accounted for, and, as the main reason, we have the isoelectronium instead 
of one single electron. 

In the Appendix, we present the results obtained from the CSF based 
recurrence relations by numerical calculations for ordinary H2 ion and for 
H2 system, at M = 2me. These results are presented in Tables ^ and ^ 
Tables ^ and ^ have been derived from Tables ^ and ^ respectively, by 
simple adding internuclear potential energy 1/R, to obtain the total energy 
of the system. 

In Table 0, we present the Isag electronic term of if 2^. In Table ^, we 
present the total energy of H2. In Table |, we present the Isag term of H2, 
at the mass M = 2me. In Table |^, we present the total energy of H2, at the 
mass M = 2me. Also, in Table ^, we present the total minimal energies of 
H2 and optimal distances R, for various values of the isoelectronium mass 
parameter, M = rjme. All the data of these Tables are purely theoretical 
and, additionally, we plot them in Figures 1-8, for the reader convenience. 
Figures 6 and 8 give more detailed view on the interval 0.26 < rj < 0.34. 

The analysis of the data in Tables |^-^ is simple. Namely, one should iden- 
tify the minimal value of the energy in each Table. One can use Figures 1-4 
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for visual identification of the minima, and then turn to the corresponding 
Tables 0-^, to reach a higher numerical accuracy. 

Let us consider, as an example. Table 13. One can see that the energy 



minimum for H2 is 



Eua = -2.0 a.u. at R = a.u. (4.1) 



Note that Eisa is the electronic energy, i.e. the internuclear repulsion has 
not been taken into account here. Our remark is that this energy value 
corresponds to He~^ ion due to the fact that two H2 nuclei are superimposed 
and form He nucleous at R = 0. 

Let us consider now Table ^. In this Table, one can find the line |2.0 | — 
0.602634| corresponding to visual minimal value of the energy. To identify a 
more precise value of the minimal energy, one should use the interpolation of 
the data. This gives us the minimum of the total energy of H2, 

E = -0.6026346 a.u. at R = 1.9971579 a.u. (4.2) 

This theoretical value represents rather accurately the known experimental 
value E'exper [-f^2^] = —0.6017 a.u. [|1^] for H2 ion, thus establishing the valid- 



ity of our CSF based calculations. For completeness, we note that the experi- 
mental dissociation energy of H2 ion is Dexper[H2] — 0.0974 a.u. = 2.65 eV, 
and the internuclear distance Rexper[H2] — 2.00 a.u. = 1.0584A. 

Let us now consider Table ^. One can see that the energy minimum for 
H2, at M = 2me, is 

Eisa = —16.0 a.u. at -R = a.u. (4.3) 

Note that Eisa only yields the isoelectronium's energy, i.e. the internuclear 
repulsion has not been taken into account here. Our remark is that this 
energy value corresponds to the He atom, where the two electrons form a 
stable point-like isoelectronium of mass M = 2m f.. 

Let us consider now Table ^ which is of striking interest for our study. 
In this Table, one can find the line |0.250 | - 7.61428940411169996] cor- 
responding to the visual minimal value of the energy. To identify a more 
precise value of the minimal energy, one should use the interpolation of the 
data given in this Table. This gives us the minimum of the total energy of 

E = -7.617041 a.u. at i? = 0.258399 a.u. (4.4) 
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This theoretical value is in quite good agreement with the preceding the- 
oretical result by Santilli and Shillady obtained via structurally different 
variational numerical method, E^ar = —7.61509174 at Ryar = 0.2592 (see 
third column of Table 1 in Ref.[|l|). It is quite naturally to observe that this 
variational energy is a bit higher (by 0.002 a.u.) than the above one obtained 
from the exact solution (as it is expected to be for any variational solution). 
However, this exact theoretical value (4^) does not meet the experimental 



value Eexper[H2] = —1.17 a.u. |jT5[ known for H2 molecule. Indeed, adopted 



approximation that the isoelectronium is point-like, stable, and has mass 
M = 2me leads us to the theoretical value ( |4.4| ) while the known experimental 
value, Eexper[H2] = —1.17 a.u., differs much from it. 

Essentially the same conclusion is due to numerical program SASLOBE 
by Santilli and Shillady |l[, where Gaussian screened Coloumb potential in- 
teraction between the electrons, rather than the stable point-like isoelec- 
tronium approximation, has been used to achieve final precise fit of = 
— 1.174474 a.u., with the obtained bond length R = 1.4011 a.u., at the iso- 
electronium correlation length Vc = 0.01125 a.u. (see Table 1 in Ref.[^). We 
discuss on this issue in Sec. |^. 

Our remark is that, due to Table |^, the experimental value E = —1.17 a.u. 
is fitted by the distance R = 0.072370 a.u. However, this energy value is 
not minimal and thus can not be ascribed reasonable physical treatment in 
Table |. 

Our conclusion from the above analysis is that we have two main possi- 
bilities to overcome this sharp discrepancy between our theoretical and the 
experimental binding energy values which has place at M = 2me: 

1. Consider unstable isoelectronium, i.e. the four-body Santilli-Shillady 
model of H2 molecule; 

2. Treat the mass M of isoelectronium as a free parameter, instead of 
fixing it to M = 2me, assuming thus some defect of mass discussed in 
Introduction, 

in order to fit the experimental data on H2 molecule. 

The first possibility will be considered in a subsequent paper because it 
needs in application of different technique, while the second possibility can 
be studied within the three-body Santilli-Shillady model under consideration 
to which we turn below. 
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In the next Section, we develop simple formalism allowing one to deal with 
the mass and charge of isoelectronium viewed as free parameters, and arrive 
at the conclusion (see Table ^) that the restricted three-body Santilli-Shillady 
model of H2 molecule is capable to fit the experimental binding energy, with 
the total mass of isoelectronium equal to M = 0.308381me, although with 
the internuclear distance about 19.6% bigger than the experimental value. 



4.2 The scaling method 



In order to relate the characteristics of H2 ion like system to that of thor- 
oughly studied H2 ion, we develop scaling method based on the Schrodinger 
equation. The neutral H2 system with stable point-like isoelectronium is an 
example of H2 ion like system in which we are particularly interested here. 
Below, we develop scaling method for the case of arbitrary mass and charge 
of the particle. 

Let us write the Schrodinger equation for a particle of the rescaled charge 

e ^ -Ce, (4.5) 
(we turn here from e = — 1 to — e = 1 representation), and the rescaled mass 

me rirrie, (4.6) 
with equal charges of nuclei, +eZi = +eZ2 = +eZ, 



(Ze' (Ze' Z^e^ 



n 



+ 



R. 



ab 



(4.7) 



where r/ and C are scaling parameters, and Rab is distance between the nuclei. 
The condition Zi = Z2 is an essential point to stress here because owing to 
which we can successfully develop the scaling method. We introduce the unit 
of length, 

1 



ro = —TT^rs = ^, (4.8) 

ri(^Z TjQZ meC^ 

where tb is Bohr's radius. Dividing Eq.( ^.7| ) by C^^^, and multiplying it by 

ro, we get 

fi^ l_o 1 1 ( Z\ 1 

Tr, n \ C J Rab 



rj^Zfriee' 



■2 ra n 



CZe^ 



ij. (4.9) 
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We introduce dimensionless entities p = r/ro, pa = Ta/ro, pi, = rb/ro, and 
TZ = Rab/ro- Then, Laplacian in Eg. ( [4 .91) becomes ^qV^ = V^. Further, 
introducing unit of energy, 

Eo 



(4.10) 



we have dimensionless energy e = E/Eq so that Eq. (^4.9|) can be rewritten as 

1 1 1 
Pa Pb (47^) 



2 P 



(4.11) 



Note that, at ?7 = 1, C = 1; and Z = 1, the constants ro{riX,Z) and 
Eq^t], (, Z) reproduce ordinary atomic units. 



ro(l, 1, 1) = rs 



rripC^ 



Eo{l, 1, 1) = 2Eb 



ruee 



h 



2 ' 



(4.12) 



and we recover the case of H2 ion. On the other hand, in terms of dimen- 
sionless entities the original Schrodinger equation for H2 ion is 



1 

pa 



1 

pb 



1 

R 



ipo = e{R)ipo, 



(4.13) 



where R = Rab/rs- Comparison of Eq.( |4.11| ) and Eq.( |4.13|) shows that by 
putting R = {(/Z)7l in Eq.( |4.11| ), we obtain the equation. 



2 ' 



1 1 1 
+ p 

Pa Pb R 



ij = e{R)ij, 



(4.14) 



which identically coincides with the original Eq. (|4.13| ). The difference is that 
Eq.( |4.14| ) is treated in terms of the rescaled units, ro(?7, C, Z) and Eo{ri, (, Z), 
instead of the ordinary Bohr's units, and Eb- As the result, we have one 
and the same form of Schrodinger equation for any H2 like system charac- 
terized by equal charges of nuclei. This makes a general ground to calculate 
some characteristic entity of any H2 like system when one knows its value 
for H2 ion. 

Particularly, one can easily derive Rab and E for the system with arbitrary 
parameters rj, (, and Z from their values, Rab[H2] and E[H2] = 2EB — e{R), 
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obtained for H2 ion (for which 7] = 1, ( = 1, and Z = 1). Indeed, since for 
arbitrary t], (, and Z 

R = ^n=^^ = ^^vCZ, (4.15) 
Z Z ro Z tb 

we can estabhsh the following relationship between the distances correspond- 
ing to arbitrary ZC^Z system and H2 ion, 

Rab = (4.16) 

It is remarkable to note that the dependence on Z disappeared in Eq. (|4.16|) . 
In the case of isoelectronium of mass M = 2me and charge — 2e, we have 
rj = 2 and C = 2, so that 

Rab = (4.17) 



Also, the energy E{R) of Z(Z system and energy e{R) of H2 ion are related 
to each other according to the equation. 



/ 4 ^ 

*2 ^2 / "^eC 



E{R)=vCz'[^UiR)- (4.18) 



4.2.1 The case M = 2me 

So, in the case of isoelectronium of mass M = 2me and charge — 2e, we get 

E{Rab) = 8eiR). (4.19) 

Note however that the factor (/Z = 2 arised due to i? = {(/Z)7l is hidden 
here so that in order to calculate the values of E{Rab) and Rab from e{R) 
and R respectively one should multiply £ by 8 and -R by 1/4. 

As the result, in accordance with the scaling method the points can be 
calculated due to the following rule: 

{R, E) ^{R, E + 1/R) (i?/4, 8E) (i?/4, 8E + 4/R), (4.20) 

for Tables |^ ^ |^ ^ ^ ^. One can easily check numerically that these prop- 
erties indeed hold true for the presented Tables. Thus, the scaling method 
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can be used instead of the independent numerical calculations for H2 system 
if one has the data for H2 ion. 

It is highly important to note here that the energy minimum in Table ^ is 
not rescaled to the energy minimum in Table ^ due to the absence of energy 
scaling between these Tables; see Eg. (H^.2CI|) , from which one can observe that 
{8E + 4/i?) can not be expressed as n{E + 1/-R), where n is a number. So 
one needs to identify minimum in Table ^ independently (after calculating 
all the points) , rather than direct rescale the minimum from Table ^ to try 
to get minimum for Table |. 



4.2.2 The case M = rjme 

For a more general case of isoelectronium mass, M = rjirie, and charge — 2e, 
we should keep the following sequence of calculations: 

(R, E) ^{R, E + 1/R) Ar^E) ^ Ar^E + |). (4.21) 

starting from Table ^ to obtain, at the last step, the table of values (similar to 
Table ^ from which we should extract a minimal value of the energy and the 
corresponding optimal distance, at each given value of mass rj. The result of 
the analysis of a big number of such tables is collected in Table ^, where the 
interval 0.26 < rj < 0.34 appears to be of interest; M = 77, in atomic units. 
Plots of the data of Table ^ are presented in Figures 5 and 7 (Figures 6 and 
8 give more detailed view on the interval of interest) show that 

E„,in{M) - -3.808M, Ropt{M) ~ (4.22) 

to a good accuracy. Note that Emin{M) unboundedly decreases with the 
increase of M (there is no local minimum), so we can use a fit, instead of the 
minimization in respect with M. From this Table, we obtain the following 
final fit of the binding energy for the restricted three-body Santilli-Shillady 
model of H2 molecule: 

M = 0.308381me, E = -1.174475 a.u., 1.675828 a.u., (4.23) 

where the mass parameter M of the isoelectronium has been varied in order to 
meet the experimental energy Eexper[H2] = —1.174474 a.u. = —31.9598 eV. 
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Using this value of mass, M = 0.308381me, we computed the total energy 
as a function of the internuclear distance R, and depicted it in Figure 9 to 
illustrate that R = 1.675828 a.u. indeed corresponds to a minimal value of 
the energy. Note that the predicted optimal distance R = 1.675828 a.u. = 
O.8868IOA appears to be about 19.6% bigger than the conventional experi- 
mental value Rexper[H2] = 1.4011 a.u. = O.742A. 

This rather big (19.6%) discrepancy can not be ascribed to the Born- 
Oppenheimer approximation used in this paper since it gives relatively small 
uncertainty in the energy value, even in the case of H2 molecule. We stress 
here that in the Born-Oppenheimer approximation, the three-body problem 
(the Schrodinger equation) can be given exact solution owing to separation 
of the electronic and nuclear degrees of freedom while the full three-body 
problem (accounting for the wave functions of the nuclei, etc.) can not be 
solved exactly. 

In a strict consideration, we should calculate the dissociation energy of H2 
molecule, D = 2Eq—E—E"'^'^\ to make comparison to the experimental value, 
Dexper[H2] ^ 0.164 a.u. = 4.45 eV Here, Eq = -0.5 a.u. = -13.606 eV 
is the ground state energy of separate if-atom and E"-"^^^ is the energy of 
zero mode harmonic oscillations of the nuclei, with the experimental value 



exper L 



E^^^L\H2\ ~ 0.01 a.u. = 0.27 eV |T5|. One can see that the zero mode 



energy E"'^'^'' (which is taken to be = 0, in the Born-Oppenheimer 

approximation) is estimated to be less than 1% of the predicted E. The 
leading anharmonic correction to the harmonic oscillation energy is estimated 
to be 1.4% of i.e. it is of the order of 0.00014 a.u. = 0.004 eV, in 

the case of H2 molecule. So, in total the Born-Oppenheimer approximation 
makes only up to 1% uncertainty, which is obviously insufficient to treat the 
predicted R = 1.675828 a.u. as an acceptable value, from the experimental 
point of view. 

Note that, at the given M, we can not "fix" R to be equal to the de- 
sired experimental value Rexper = 1.4011 a.u. unless we shift E to some 
nonminimal value, which is, as such, meaningless. Conversely, if we would fit 
experimental Rexper by varying M, we were obtain Emin between —1.52 a.u. 
and —1.33 a.u. (see Table |^), which is much deviated from the experimen- 
tal Eexper- In othcr words, the relation between E and R, governed by the 
Schrodinger equation, is such that at some value of R there is a minimum of 
E so that R is not some kind of free parameter here since the system tends 
to minimize its own energy. In accordance to the exact solution of the model. 
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our single free parameter, M, can not provide us with the exact fit of both 
the experimental values, E^xper and Rexper- 

Thus, we arrive at the conclusion that the three-body Santilli-Shillady 
model of H2 molecule yields the result ( ^4.23|) , which indicates that the as- 



sumption of stable point-like isoelectronium builds a crude approximation 
to the general (four-body) Santilli-Shillady model. This means that we are 
forced to possess that the isoelectronium is not stable point-like quasi-particle, 
to meet the experimental data on H2 molecule. 

4.3 Variational solution 

In studying H2 ion like systems, one can use Ritz variational approach to 
obtain the value of the ground state energy as well. This approach assumes 
analytical calculations, which are easier than that used in finding the above 
exact solution but they give approximate value of the energy. It is helpful in 
making simplified analysis of the system. This can be made for the general 
case of isoelectronium total mass, which eventually undergoes some " defect" 
while its "bare" total mass is assumed to be M = 2me. Ritz variational solu- 
tion of the H2 like problem yields, of course, similar result for the energy of 
H2. Below, we present shortly results of our calculations. However, we stress 
that the variational solution is given here just to make some support to the 
exact solution, and to see the typical order of the variational approximation. 

Using hydrogen ground state wave function and one-parameter Ritz vari- 
ation, we obtain the following expression for the energy of H2 like system: 

E(p) = + e^l Jl + p)e-^^ + il-y)e-^ 

2ao aop 1 + (1 + p + ip2)e-P ' ^ ' ^ 

where p = R/oq is variational parameter. For the general case of mass m 
and charge q = (e, Eg. ( [4. 241) can be rewritten in the following form: 



E{p,0 = ^^[-^ + F{p)), (4.25) 

where 

^^'^ - -p l + (l + p+ip2)e-P ■ ^'-2^) 
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Numerically, the function F{p) reaches minimum at the value p = 2.5, which 
should be used in the above expression for E{pX)- So, putting ^ = 1 we 
obtain the variational value of H2 ion energy, E{p) = —0.565. Note, to 
make a comparison, that we have the value E^xact = —0.6026 due to the 
exact solution (|4.2| ), and the value Eexper[H2] = —0.6017 as the experimental 
value of the energy of H2 ion. Thus, the optimal distance between the 
protons in H2 ion is Rm = clqP = 2.5 a.u., and the obtained variational 
energy E is slightly higher than both the values E^xact and Eexper, as it is 
normally expected to be in the variational approach. Now, we should replace 
electron by isoelectronium to describe the associated H2 model. Substituting 
M = 2me and C = 2, we see that the r.h.s. of Eg. ( [4. 251) contains overall factor 
8, in comparison to the H2 ion case (M = rUe and ( = 1), 

E{p)=8\2EB\(-^ + Fip)y (4.27) 

The function F{p) remains the same, and its minimum is reached again at p = 
2.5. Then, energy of H2 molecule due to Eq.( [4.27| ) is E{p) = — 8|2i?B|0.565 = 
—4.520 a.u. This value should be compared with the one given by Eg. (|4.4]) . 

Below, we collect the above mentioned data and results of this Section in 
Table |. 





E, a.u. 


R, a.u. 


H2 ion, exact theory (A^=16) 
H2 ion, experiment ||l5tl 


-0.6026346 
-0.6017 


1.9971579 
2.00 


3-body H2, M=2me, exact theory (A^=16) 
3-body H2, M=2me, var. theory Q 

3- body H2, M=0.381me, exact theory (iV=16) 

4- body H2, rc=0. 01125 a.u., var. theory ^ 
H2, experiment 


-7.617041 

-7.61509174 

-1.174475 

-1.174474 

-1.174474 


0.258399 

0.2592 

1.675828 

1.4011 

1.4011 



Table 1: The total ground state energy E and the internuclear distance R. 
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5 Concluding remarks 



In this paper we have shown that the restricted three-body SantiUi-Shillady 

isochemical model of the hydrogen molecule admits an exact analytic solu- 
tion capable of representing the molecular binding energy in a way accu- 
rate to the sixth digit, E = —1.174475 a.u., and the internuclear distance 
R — 1.675828 a.u., which is about 19.6% bigger than the conventional ex- 
perimental value, Rexper[H2\ — 1.4011 a.u. 

We should emphasize that the presented exact analytical solution includes 
infinite chain fractions. They still need numerical computation to reach the 
characteristic values of H2 ion like systems, such as the ground state energy, 
with the understanding that these values can be reached with any needed 
accuracy. For example, at the lengths of the chain fractions = 100 and 
N — 50 for the angular and radial eigenvalues, one achieves accuracy of the 
ground state energy of about 10^^^. 

The general (four-body) SantiUi-Shillady isochemical model of H2 cannot 
be, apparently, solved exactly, even in Born-Oppenheimer approximation, so 
that Ritz variational approach can be applied here to get the approximate 
values of the ground state energy and corresponding internuclear distance. 

Ritz variational approach is a good instrument to analyze few-body prob- 
lems, and restricted H2 molecule is such a system. It is wellknown that the 
variational solution of the ordinary model of H2 molecule includes rather com- 
plicated analytical calculations of the molecular integrals, with the hardest 
part of work being related to the exchange integral. Particularly, evaluated 
exchange integral for H2 molecule is expressed in terms of a special func- 
tion (Sugiura's result, 1927). It is quite natural to expect that even more 
complications will arise when dealing with the Hulten potential. 

The reason to consider the general four-body SantiUi-Shillady model of H2 
molecule, after the made analysis of H2 like system approximate approach 
to it, is that the stable point-like isoelectronium H2 model-based theoreti- 
cal prediction does not meet the experimental data on H2 molecule, for the 
"bare" isoelectronium mass M — 2me, although we achieved essentially ex- 
act representation of the binding energy taking M — 0.308381me. Also, this 
stable point-like isoelectronium (three-body) model does not account for es- 
sential effect existing in the general (four-body) model. This effect is related 
to the potential barrier between the region associated to the attractive Hulten 
potential, ri2 < tq, and the region associated to repulsive Coloumb poten- 
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tial, ri2 > ro, where tq is the distance between the electrons at which Hulten 
potential is equal to Coloumb potential, V{ri2) = 0; see Eq.( p.5|) . Char- 
acteristics of the barrier can be extracted from the function V{ri2). The 
barrier is finite for the used values of the parameters Vq and Vc so that the 
electrons penetrate it. The two Is electrons are thus simultaneously in two 
regimes, the first is strongly correlated regime due to short-range attractive 
Hulten potential (isoelectronium) and the second is weakly correlated regime 
due to the ordinary Coloumb repulsion. Also, there exist a transient regime 
corresponding to the region about the equilibrium point, ~ ro, i.e. inside 
the barrier. Schematically, one could thought of that the electrons are, for 
instance, 10% in the isoelectronium regime, 1% in the transient regime, and 
89% in the Coloumb regime. We stress that in the three-body approach to 
H2 molecule considered in this paper we have 100% for the isoelectronium 
regime. 

Numerical computation by Santilli and Shillady [|1| based on Gaussian 
transform techniques and SASLOBE computer program has shown excel- 
lent agreement of the general four-body model with experimental data on 
H2 molecule. They used Gaussian screened Coloumb potential as an ap- 
proximation to the Hulten potential. It would be instructive to use Ritz 
variational approach, which deals with analytical calculations, in studying 
the four-body Santilli-Shillady model of H2 molecule. One can try it first for 
the Gaussian screened Coloumb potential, or exponential screened Coloumb 
potential in which case there is a hope to achieve exact analytical evalua- 
tion of the Coloumb and exchange integrals. Being a different approach, this 
would give a strong support to the numerical results on the ground state 
energy obtained by Santilli and Shillady. Also, having analytical set up one 
can make qualitative analysis of the four-body Santilli-Shillady model of H2 
molecule. However, we should to note that these potentials, being approxi- 
mations to the Hulten potential, will yield some approximate models, with 
corresponding approximate character of the results. 
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Appendix 



We use = 16 power degree approximation, the polynomials Qtv^ and 

to find both the radial, X^^\p,a), and angular, \^'^\p,b), eigenvalues of the 

CSF. Qn^s are obtained by the use of recurrence relations ( p.76| ) and defini- 



tions of the coefficients as, Ps, Is, ps, ^s, and 5s, where we put 5 = 0, i.e. 
Zi = Z2 = 1, and quantum number m = 0. Each of the two polynomials 
has 16 roots for A from which we select one root which is appropriate due 
to its asymptotic behavior at i? 0. Numerical solution of the equation 
X^^\p, a) = \^y\p, h) gives us the list of values of the electronic ground state 
energy E{R) = Eisa{R), which corresponds to Isag term of the H2 ion, as 
a function of the distance R between the nuclei. Table ^ presents the result, 
where no interpolation has been used. Numerical computation of each point 
in Table |^ took about 88 sec on ordinary Pentium desktop computer. 

Below, we present some useful numerical values enabling one to convert 
atomic units, at which mg = e = ^ = 1, to the other units. Note also that 
for the energy 1 a.u. = 1 hartree, and for the length 1 a.u. = 1 bohr. 

Atomic units in terms of the other units 



1 a.u. of mass, rrie 


9.10953-10-^« 


gramms 


1 a.u. of charge, e 


1.60219-10-1^ 


Coloumbs 


1 a.u. of action, h 


1.05459-10-27 


erg - sec 


1 a.u. of length, 


0.529177-10-^ 


cm 


1 a.u. of energy, 


27.2116 


eV 


1 a.u. of time, 


2.41888-10-1^ 


sec 


1 a.u. of velocity, e'^/h 


2.18769-10^ 


cm/sec 


a = f- 

nc 


1/137.0388 





Conversion of the energy units 





a.u. 


eV 


Kcal-mole 


cm 1 


a.u. 


1 


27.212 


6.2651-102 


2.1947-10^ 


eV 


3.6749-10-2 


1 


23.061 


8065.48 


Kcal-mole 


1.5936-10-3 


4.336410-2 


1 


3.4999-102 


cm"! 


4.5563 -10-6 


1.239810"^ 


2.8573-10-3 


1 
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R, a.u. E{R), a.u. 



R, a.u. E{R), a.u. 



0.0 


-1 


99999999761099225 


2.01 


-1 


10013870349441877 


0.1 


-1 


97824134920757046 


2.05 


-1 


09030214496519061 


0.2 


-1 


92862028526774320 


2.1 


-1 


07832542203506842 


0.3 


-1 


86670393395684293 


2.2 


-1 


05538508113994433 


0.4 


-1 


8UU754U5253452878 


2.3 


-1 


03371349485820318 


0.5 


-1 


73498799160719041 


2.4 


-1 


01322030525887973 


0.6 


-1 


67148471440012302 


2.5 


-0 


99382351101203490 


0.7 


-1 


61119720301672586 


2.6 


-0 


97544858094023219 


0.8 


-1 


55448006449772595 


2.7 


-0 


95802766004904907 


0.9 


-1 


50138158334624467 


2.8 


-0 


94149886061738322 


1.0 


-1 


45178630031448951 


2.9 


-0 


92580563147803989 


1.1 


-1 


40550252191841256 


3.0 


-0 


91089619738235434 


1.2 


-1 


36230785783351171 


3.1 


-0 


89672306127076382 


1.3 


-1 


32197139010318509 


3.2 


-0 


88324255989570446 


1.4 


-1 


28426925496894185 


3.3 


-0 


87041447461742614 


1.5 


-1 


24898987186705512 


3.4 


-0 


85820167794222435 


1.6 


-1 


21593722446146546 


3.5 


-0 


84656982450508629 


1.7 


-1 


18493139974611416 


3.6 


-0 


83548707392472181 


1.8 


-1 


15580915764590441 


3.7 


-0 


82492384412924800 


1.9 


-1 


12842156954614458 


3.8 


-0 


81485259165546253 


1.95 


-1 


11533575206408963 


3.9 


-0 


80524772550750985 


1.99 


-1 


10514450160298682 


4.0 


-0 


79608496995054425 


2.0 


-1 


10263415348745197 


8.0 


-0 


62757022044109352 



Table 2: The electronic energy of H2 ion (see Fig. 1). 



Minimum of the energy EuaiR) is E^a = -1.9999999976 a.u. at i? = 0, 
which reproduces the known value Eisa = —2 a.u. to a very high accuracy. 
Moreover, one can compare Table |^ and the table of Ref . to see that each 
energy value in Table |^ does reproduce Wind's result up to seven decimal 
places. This means that our numerical calculations are correct. 

We remark that Wind used = 50 approximation and presented seven 
decimal places while we use N = IQ approximation and present seventeen 
decimal places. Alas, there is no need to keep such a high accuracy, and also 
Wind mentioned that even = 10 approximation gives the same result, up 
to seven digits. 

By making 16th-order interpolation of the points in Table |^ and adding 
to it the potential of interaction between the nuclei, 1/R,we obtain the list of 
values of the total energy presented in Table |^. It reveals the only minimum 
of the total energy, E{R) + R~^ = Emin = —0.6026346 a.u. at the distance 
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R = R , = 1.9971579 a.u. 



Pan 

-ft, a.u. 


-Cjyrt) -\- rL , a.u. 


P Q 11 

it, a.u. 


ii/\rLj -\- it , a.u. 


0.1 


A-R 091 7fi 


2.1 


U.UU^J-Ot-l 


0.2 


+3.07138 


2.2 


-0.600840 


0.3 


+1.46663 


2.3 


-0.598931 


0.4 


+0.69924 


2.4 


-0.596554 


0.5 


+0.26501 


2.5 


-0.593824 


0.6 


-0.004818 


2.6 


-0.590833 


0.7 


-0.182626 


2.7 


-0.587657 


0.8 


-0.304480 


2.8 


-0.584356 


0.9 


-0.390270 


2.9 


-0.580978 


1.0 


-0.451786 


3.0 


-0.577563 


1.1 


-0.496412 


3.1 


-0.574142 


1.2 


-0.528975 


3.2 


-0.570743 


1.3 


-0.552741 


3.3 


-0.567384 


1.4 


-0.569984 


3.4 


-0.564084 


1.5 


-0.582323 


3.5 


-0.560856 


1.6 


-0.590937 


3.6 


-0.557709 


1.7 


-0.596696 


3.7 


-0.554654 


1.8 


-0.600254 


3.8 


-0.551695 


1.9 


-0.602106 


3.9 


-0.548837 


2.0 


-0.602634 


4.0 


-0.546085 



Table 3: The total energy of H2 ion (see Fig. 2). 



The results collected in Table ^ have been obtained directly by numerical 
calculations with the use of replacements and a 4a, where p 

and a are defined by Eg. (|3.5CI| ), in the coefficients a^, [3s, 7s, Ps, and 
5s of the recurrence relations. These replacements have been made due to 
Eq. (|3.15|) , with the mass parameter M = 2 and the charge parameter q = —2, 
corresponding to the stable point-like isoelectronium of mass M = 2me and 
charge — 2e. In addition, it turns out that Table ^ can be derived directly from 
Table |^ by the use of rescalements R — > R/4 and E 8E. This remarkable 
property is confirmed by the scaling method developed in Sec. [4.2| , and proves 
that the scaling method is correct. By adding 1/Rto the isoelectronic energy 
values of Table ^ we obtain Table ^ showing the total energy of the H2 
system, at the mass M = 2me. The minimum of the total energy is found 
E{R) + R-^ = Emin = -7.617041 a.u. at R = Ropt = 0.258399 a.u. 
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R, a.u. E{R), a.u. | R, a.u. E{R), a.u. 






00 


-16.0000000000000008 





50 


-8 


82107371596060652 





05 


-15.4289613540288446 





55 


-8 


44308064942978475 





10 


-14.4060324481253427 





60 


-8 


10576243870545276 





15 


-13.3718774940028106 





65 


-7 


80358864752256309 





20 


-12.4358407555134897 





70 


-7 


53199088488360768 





225 


-12.0110526341241952 





75 


-7 


28716957910914597 





25 


-11.6142894041116995 





80 


-7 


06594047783391143 





275 


-11.2440221722385613 





85 


-6 


86561342351093717 





30 


-10.8984629511501962 





90 


-6 


68389659171409977 





35 


-10.2741538838779300 





95 


-6 


51882073404630535 





40 


-9.72749779604447084 


1 


00 


-6 


36867910416260141 





45 


-9.24647351730713928 











Table 4: The isoelectronium energy of the three-body H2 system, at the mass 
M = 2me (see Fig. 3). 



R, a.u. E{R), a.u. | R, a.u. E{R), a.u. 






010 


-1-84.0273769895390465 





150 


-6 


70521082733614370 





015 


-1-50.7305619003569852 





200 


-7 


43584075551348888 





020 


-1-34.1124961901906909 





250 


-7 


61428940411169996 





025 


-1-24.1713076006601701 





300 


-7 


56512961781686321 





030 


-t-17.5720885190207170 





350 


-7 


41701102673507239 





035 


-H2.8849357295675638 





400 


-7 


22749779604447084 





040 


-1-9.39442987714794597 





450 


-7 


02425129508491696 





045 


-1-6.70277469374505585 





500 


-6 


82107371596060652 





050 


-1-4.57103864597115538 





550 


-6 


62489883124796641 





055 


-1-2.84698239022676524 





600 


-6 


43909577203878669 





060 


-1-1.42896579635860376 





650 


-6 


26512710906102388 





065 


-1-0.24650856005471055 





700 


-6 


10341945631217797 





070 


-0.75082137484248256 





750 


-5 


95383624577581205 





075 


-1.60007720159364552 





800 


-5 


81594047783391143 





080 


-2.32910343012903808 





850 


-5 


68914283527564279 





085 


-2.95923562411244667 





900 


-5 


57278548060298906 





090 


-3.50710131449860362 





950 


-5 


46618915509893721 





095 


-3.98585308299203866 


1 


000 


-5 


36867910416260141 





100 


-4.40603244812534278 











Table 5: The total energy of the three-body H2 system, at the mass M = 2me 
(see Fig. 4). 
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IVl , a.u. 


^mml^Wj, a.u. 


R-opti-'y-l ) , a.u. 


0.10 


-0.380852 


5.167928 


0.15 


-0.571278 


3.445291 


0.20 


-0.761704 


2.583964 


0.25 


-0.952130 


2.067171 


0.26 


-0.990215 


1.987664 


0.27 


-1.028300 


1.914050 


0.28 


-1.066385 


1.845688 


0.29 


-1.104470 


1.782044 


0.30 


-1.142556 


1.722645 


0.307 


-1.169215 


1.683367 


0.308 


-1.173024 


1.677899 


0.308381 


-1.174475 


1.675828 


0.309 


-1.176832 


1.672471 


0.31 


-1.180641 


1.667073 


0.32 


-1.218726 


1.614977 


0.33 


-1.256811 


1.566041 


0.34 


-1.294896 


1.519981 


0.35 


-1.332982 


1.476553 


0.40 


-1.523408 


1.291982 


0.45 


-1.713834 


1.148428 


0.50 


-1.904260 


1.033585 


0.75 


-2.856390 


0.689058 


1.00 


-3.808520 


0.516792 


1.25 


-4.760650 


0.413434 


1.50 


-5.712780 


0.344529 


1.75 


-6.664910 


0.295310 


2.00 


-7.617040 


0.258396 







Table 6: The minimal total energy i?mm and the optimal internuclear distance 
Ropt of the three-body H2 system as functions of the mass M of the stable 
point-like isoelectronium (see Figs. 5-8). 



Table ^ presents result of calculations of the minimal total energies and 
corresponding optimal distances, at various values of the isoelectronium mass 
parameter M = rjirtf. (M = rj, in atomic units). We have derived some 27 
tables (such as Table ^ from Table |] by the scaling method according to 
Eq.( |4.21| ), and find minimum of the total energy in each table, together 
with the corresponding optimal distance. Then we collected all the obtained 
energy minima and optimal distances in Table ^. With the fourth order in- 
terpolation/extrapolation, the graphical representations of Table ^jshow (see 
Figures 5-8) that the minimal total energy behaves as Emin{M) ~ — 3.808M, 
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and the optimal distance behaves as Ropt{M) ~ 0.517/M, to a good accu- 
racy. One can see that at M = 2me we have Emin{M) = —7.617040 a.u. 
and Ropt{M) = 0.258396 a.u., which recover the earher obtained values 
Emin = —7.617041 a.u. and Ropt = 0.258399 a.u. of Table ||, to a high ac- 
curacy, thus showing once again correctness of the used scaling method. In 
fact, the values of E and R for M = 1.50me, M = 1.75me, and M = 2.00me 
in Table ^ have been obtained by extrapolation so they are not as much ac- 
curate as they are in Table ^. However, this is not of much importance here 
because we use them only to check the results of the scaling method. 

The main conclusion following from Table |^ is that the mass parameter 
value M = 0.308381me fits the energy value E^in{M) = -1.174475 a.u., 
with the corresponding Ropt{M) = 1.675828 a.u., which appears to be about 
19.6% bigger than the experimental value Rexper[H2] = 1.4011 a.u. The total 
energy as a function of internuclear distance, for this value of mass, M = 
0.308381me, is shown in Figure 9 to illustrate that the obtained optimal 
distance Ropt = 1.675828 a.u. corresponds to a minimal value of the total 
energy. 
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12 3 4 

R, a . u . 




Figure 2: The total energy E{R) + R'^ of Ht ion as a function of the 
internuclear distance R. 
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Figure 3: The isoelectronium energy E{R) of the H2 system as a function of 
the internuclear distance R, at the isoelectronium mass M — 2me. 




0.2 0.4 0.6 0.8 1 

R, a.u. 



Figure 4: The total energy E{R) + R ^ oi the H2 system as a function of the 
internuclear distance R, at the isoelectronium mass M — 2me. 
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0.5 1 

Mass M, a.u. 



Figure 5: The minimal total energy £^^j„(M) of the H2 system as a function 
of the isoelectronium mass M. 




-1.3 



0.26 0.28 0.3 0.32 

Mass M, a.u. 



0.34 



Figure 6: The minimal total energy £'^j„(M) of the H2 system as a function 
of the isoelectronium mass M. More detailed view. 
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5^ \ 




0.5 1 1.5 2 

Mass M, a.u. 



Figure 7: The optimal internuclear distance Ropt{M) of the H2 system as a 
function of the isoelectronium mass M. 




0.26 0.28 0.3 0.32 0.34 

Mass M, a.u. 



Figure 8: The optimal internuclear distance Ropt{M) of the H2 system as a 
function of the isoelectronium mass M. More detailed view. 
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